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w^ . Abstract 

. ~^ ' In recent years, several particle-based stochastic simulation algorithms (PSSA) 

i-Q . have been developed to study the spatially resolved dynamics of biochemical net- 

^H' works at a molecular scale. A challenge all these approaches have to address is 

to allow for simulations at cell-biologically relevant timescales without neither 
neglecting important spatial and biochemical properties of the simulated system 
►^ . nor introducing ad-hoc assumptions not based on physical principles. Here we 

f^ ' describe a PSSA that permits large time steps while still retaining a high de- 

C^ . gree of accuracy. The approach addresses the typical disadvantage of Brownian 

Cn ' dynamics, namely the need to use small time steps to resolve bimolecular en- 

counters accurately, by estimating the number of otherwise unnoticed encounters 



q 

t^^ ' with the help of the Green's functions of the diffusion equation incorporating 

^— ^ . molecular interactions. This method has previously been proposed for purely 

absorbing boundary conditions and irreversible bimolecular reactions. Building 
on those ideas, we developed a general-purpose PSSA that is applicable to a 
broad class of reaction-diffusion problems by incorporating reflective and radia- 
^^ , tion boundary conditions and reversible reactions. We furthermore discuss how 

Vh ' reaction-diffusion systems on 2D membranes can be described and derive small 

^ time expansions of the Green's functions that substantially speed up key cal- 

culations, particularly in the problematic case of molecules in close proximity. 
Finally, we point out the formal relationship between our and exact algorithms. 
The proposed algorithm may serve as an easily implementable and flexible, com- 
putationally efficient, coarse-grained description of reaction-diffusion systems in 
2D and 3D that nevertheless provides a stochastic, detailed representation at the 
level of individual particle trajectories in space and time. 
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1 Introduction 

Stochastic fluctuations are inherent to any ceUular biochemical system due to its 
molecular constituents. Under certain conditions, which often include the exis- 
tence of sufficiently separated timescales between different levels of organization 
of cellular biochemistry, the fluctuations can be averaged out and one arrives 
at an effective theory that treats the system at a coarser level of description 
at which fluctuations are negligible. The relevant level of description is defined 
by the experimental context. However, if the effect of fiuctuations can propa- 
gate through the different scales of the system, or if the experimental context 
changes, a theoretical treatment that explicitly takes into account fluctuations 
becomes necessary. The chemical master equation (CME) [5T], which abandons 
the notion of molecule concentrations and incorporates stochastic fluctuations, 
provides the appropriate theoretical framework. Unfortunately, it is difficult if 
not impossible to solve the CME for all but the simplest reaction networks. The 
stochastic simulation algorithm (SSA) or " Gillespie algorithm" [13 [H] avoids 
the need of finding a solution by sampling numerical realizations of the network 
components' time evolution according to the CME . However, intracellular bio- 
chemical networks operate in an environment that exhibits a spatially intricate, 
highly heterogeneous organization with aspects such as compartmentalization 
and scaffolding playing an important role. Thus, simulation tools are needed 
that are capable of abandoning the requirement of a well-stirred and homoge- 
neous environment. The SSA incorporates fluctuations but still assumes that 
the system is well-stirred and spatially homogenous. These shortcomings can be 
overcome by the so-called spatial Gillespie algorithms which partition the reac- 
tion volume in compartments small enough that within them the assumption of 
a well-stirred and homogeneous system are justified again [3J[5ni[Mj- However, 
this approach relies on the existence of a length and time scale on which the 
system is again homogeneous. Such a scale may not always exist for a particular 
intracellular reaction-diffusion network. Particle-based methods are capable of 
taking into account fluctuations and spatial aspects without introducing ad-hoc 
spatial discretizations. Such methods treat biochemical networks as composed of 
elementary unimolecular {A — >■ products) and bimolecular (A + B ^ products) 
reactions. A bimolecular A + B ^ C reaction may be depicted as a two-step 
process [351132] : 

A + b'^^A-.B^C (1.1) 

According to this kinetic scheme, a prerequisite for the occurrence of a bimolec- 
ular reaction is that the two molecules encounter each other through diffusive 
(Brownian) motion and form an encounter complex A : B. The actual chemical 
reaction is characterized by the rate constant ki. Alternatively, the encounter 
complex may decay and the molecules escape from each other with the diffu- 
sional backward rate fc_. Finally, the reaction product C can dissociate to the 
encounter complex with the dissociation constant fc_i . If the time scale associated 
with diffusion is comparable or larger than the time scale associated with the 
intrinsic chemical reaction, the bimolecular reactions are referred to as diffusion- 
influenced or diffusion-limited [32j . In these cases the Brownian motion of the 
individual molecules becomes an important element of the theoretical descrip- 
tion. Conceptually, there are two different approaches to describing Brownian 



dynamics (BD). The first approach is based on stochastic differential equations, 
often referred to as Langevin equations providing a time stepping procedure 
to generate particular stochastic realizations of a particle's space time trajec- 
tory. The second approach, in contrast, uses the associated Fokker-Planck (FP) 
equation 33 that governs the deterministic time evolution of the conditional 
probability density function (pdf) p(r, t|ro,io)- The major advantages of the 
Langevin description are its versatility and the ease of implementing it as a 
stochastic simulation algorithm. However, the Langevin method suffers from a 
severe drawback: Small time steps are in general necessary to resolve the en- 
counter events with sufficient accuracy |!6 . The need for tiny steps becomes 
especially painful in dilute systems. In these systems, the vast majority of time 
updates perform purely diffusional steps that do not involve reactions. In prin- 
ciple, this problem can be addressed by the second approach based on the FP 
equation . The Green's functions of the FP equation contain all relevant in- 
formation about the time evolution of the distribution of the particles in the 
system and the FP equation approach has been widely used in the analysis of 
diffusion- influenced reactions, cp., for instance, 0|3]. Unfortunately, for many- 
particle systems analytical solutions are not available. In addition, for more 
complicated systems it is not feasible either to generate a solution of the Fokker- 
Planck equation by numerical methods for treating partial differential equations. 
In contrast, in such systems the associated stochastic differential equation can 
still be integrated numerically. Indeed, following the seminal work by Ermak and 
McCammon [131 many simulation algorithms based on BD have been employed 
in a variety of fields, including the study of protein-protein association reactions 
[13] , dynamics of polymeric fluids [3S] and more recently, the simulation of net- 
works of interacting biomolecules [51 |3S1 123] . 

In recent years Green's function based methods have been developed that avoid 
the problems caused by many-particle systems by determining the time step 
for every simulation cycle based on the requirement that within At at most 
two particles may encounter each other, thereby excluding N-body interactions 
|38[ 139) . Factorizing a many-particle system into several independent 1-body 
and 2-body systems that can be propagated according to the Green's functions 
of the 1-body and 2-body diffusion equation permits performing large simulation 
time steps when the molecular concentrations are small. In a similar spirit, the 
first-passage kinetic Monte Carlo [IHl HZ] reduces the many-particle system to 
a set of 1-body and 2-body problems by partitioning the reaction volume into 
protective domains which contain at most two particles. Within the domains the 
particles can be propagated by first-passage and no-passage Green's functions. 
Both methods are more efficient for dilute system than conventional BD sim- 
ulations. Their event-driven nature ensures that every encounter is taken into 
account even when using large time steps. While this adds to the accuracy of 
the algorithms, for higher concentrations and particles close to boundaries the 
time step can become very small, rendering the methods inefficient. 
Here, we propose a stochastic simulation framework for reaction-diffusion pro- 
cesses that combines elements of both, the Langevin and Fokker-Planck method. 
Briefly, the particles' displacements are sampled according to the free-space 
(overdamped limit) of the Langevin equation 

^ = ;^(F + FBrownW). (1.2) 



where D,kB,T are the diffusion constant, the Bohzniann constant and the abso- 
lute temperature, respectively. The stochastic force contribution is characterized 

by 

(FBrow„)=0 and {FBrownAto)FBrown.jit))^2r)^DS,,S{to~t), (1.3) 

where ij = -jy- and the deterministic force contribution is described by F. To 
correct for underestimating the number of encounters we employ the 3D and 
2D analytical representations of the fundamental solutions of associated Smolu- 
chowski equation with absorbing, reflecting and radiation boundary conditions 
(bcs). These Green's functions can be used to compute the probability that 
particles that do not overlap after a time step nevertheless reacted at some in- 
termediate time during the time step. We emphasize that the general idea to 
use Green's functions of Fokker-Planck equations to enhance BD simulations is 
far from new, cp. [23l|2a|25l|26l[18l|30l[IIl[29]. Indeed, the method has been 
proposed for the case of purely absorbing boundary conditions and irreversible 
reactions already in [5] . Here, we build on this approach to develop a general pur- 
pose simulation framework applicable to a broad class of reaction-diffusion sys- 
tems by extending it to include reflective and radiation boundary conditions and 
reversible reactions. Furthermore, we describe how to treat reaction-diffusion 
systems in 2D. In this case, the radial Green's functions cannot be expressed by 
elementary functions, in contrast to their 3D counterparts. As a consequence, 
their numerical approximation is more costly, an issue that becomes worse for 
smaller time steps. To address this problem, we derive small time expansions for 
the key expressions. These small time expansions permit to circumvent a numer- 
ical integration and should prove useful for any simulation algorithm employing 
2D Green's functions. Finally, we establish a connection between the presented 
algorithm and the exact first-passage time algorithm to elucidate in what sense 
it can be understood as coarse-grained version. In this context an anology to the 
Gillespie formalism is helpful: While the original Gillespie algorithm is event- 
driven and takes every reaction event into account, resulting in possibly tiny 
time steps, and therefore rendering it too slow for networks with high molecular 
abundances, the tau-leaping method [T7] works with a larger (constant) time 
step At, but it has to provide a procedure to estimate the number of events that 
occur during Ai. 

Since the time step that can be taken with our approach depends less sensi- 
tively on the concentrations of the chemical species, the method is also applica- 
ble to systems for which other stochastic methods, for instance GFRD |38[I5U] . 
would fail. Compared to those approaches, the time step can be taken much 
larger while maintaining a high degree of accuracy 5:. Furthermore, chemical 
reactions near a reflecting boundary (such as a membrane) pose no particular 
challenge, again in contrast to the mentioned event-driven methods. Due to its 
well-deflned relationship with event-driven simulation algorithms, cp. section [5l 
our approach may serve as a complement for situations in which those algorithms 
are less efficient, for instance when high local particle densities would otherwise 
require using BD simulation steps 37 . 

A note on terminology: In the following we will always consider the over- 
damped limit of the Langevin equation, i.e. only position variables, but no 
velocity variables are taken into account. The corresponding FP equation is 
referred to as Smoluchowski equation. Furthermore, we will also neglect deter- 



ministic force contributfons. In this case the Smoluchowski equatfon takes the 
form of Einstein's diffusion equation [12 and in the folfowing we will use the 
terms diffusion, Smoluchowski and Fokker Planck equation interchangeably. We 
emphasize that this does not mean a limitation of the approach. Rather, as 
pointed out by [B], treating the stochastic and deterministic force contributions 
separately, offers the advantage that it is sufficent to deal only with the Green's 
functions of the diffusion equation. Thus, one can avoid using the Green's func- 
tions that take into account the deterministic force and for which only rarely 
analytical representations can be obtained. We will return to this point in sec- 
tion [lX5l 

2 Theory 

2.1 Smoluchowski equation 

The Smoluchowski equation describes Brownian motion in terms of probability 
density functions p(r, i|ro, tg). The proposed algorithm is based on the possibility 
to describe an isolated pair of diffusing particles that may react with each other 
upon encounter as the diffusion of a point-like particle near a boundary [3l [6l 
[55] . More precisely, one considers two spherical diffusing molecules A and B. 
Their probability density function is described by the two-body Smoluchowski 
equation: 

d 

—p{rA,rB,t\rAo, vbo, to) = {Da'^a + DB'^%)p{rA,rB,t\rAo, rso, ^o) (2.1) 

where Da, Db denote the diffusion constants of molecule A and B, respectively. 
By transition to the coordinates 

R - ^^A + ^VB (2.2) 

r = tb-ta (2.3) 



we obtain 



— p(R, r, i|R, r, to) - (DnVl + i?effV?)p(R, r, i|R, r, to), (2.4) 



where 



Deff = Da+Db (2.5) 

DaDb ,„ „^ 

Dr == —^ ■ (2.6) 

Evidently, the two-body Smoluchowski equation (|2.ip governs two independent 
random processes. Hence, the above equation for p(R, r, i|R, r, to) may be rewrit- 
ten as one equation for R and one for the relative inter-particle vector r 

— p(R,t|Ro,io) = i?RV^p(R,i|Ro,io) (2.7) 

d 

—p{r,t\ro,to) = Deff\/lp{r,t\ro,to), r > Oeff. (2.8) 



Ooff denotes the encounter radius and - in the absence of long range interaction 
potentials - would be given by the sum of the molecules' radii. Note that equation 
(|2.8I) describes a single particle diffusing with D^e that is excluded from a sphere 
with radius Ccff located at the origin. We will return to this analogy in section 

Defining the probability flux by 

j(r, i|ro, io) := -i?cffVrP(r, t\ro, to) (2.9) 

the FP equation (j2.8p may be written as continuity equation: 

— p(r,i|ro,to)+Vr-j(r,t|ro,io)=0, r > Oeff. (2.f0) 

The FP equations have to be completed by specifying boundary conditions for 
the conditional pdf p{r, i|ro, io) and/or the probability flux (|2.9p . Together with 
the following initial 

p(R,io|Ro,io) = <5(R-Ro) (2.11) 

and boundary condition 

p(|R| -)■ oo, i|Ro, to) = (2.12) 

equation (j2.7p is equivalent to the free-space diffusion equation with the familiar 
solution 

p(R,t|Ro,io) = ^,^nj.,,)r^. e-^^^^, (2-13) 

also known as the free-space Green's function. Henceforth, Green's functions 
will be denoted as G'(r,i|ro,io). 

The equation for the inter-particle vector r is only defined for r > acs and 
one has to impose a boundary condition for |r| = acs specifying the physics at 
the encounter distance. We will discuss the following cases ^: 

• Absorbing boundary conditions: The molecules react instantaneously upon 
encounter. 

G(|r|=aeff,t|ro,io) = 0; (2.14) 

• Reflective boundary conditions: The molecules never react upon encounter. 

J(|r|=aeff,i|ro,io)-0. (2.15) 

Here we defined the total flux through a sphere by J(r, t\ro, to) — LOd r'^^^n(r)- 
j(r, t\ro, to), d denotes the number of space dimensions, uJd '■= 2iT'^^^/T{d/2) 
is the surface area of the unit sphere in d dimensions and n(r) denotes the 
unit outward directed normal vector at a boundary point |r| = aog. 

• Radiation boundary condition [9]: Upon encounter the molecules undergo 
a chemical reaction with a certain probability, otherwise they get refiected. 

J(|r| = fleff, i|ro, to) = KaG{\r\ = Oeff, t\ro,to). (2.16) 

The constant parameter Kq relating the flux of probability to the proba- 
bility that the particles are in contact is referred to as intrinsic association 
rate at contact. It is distinct from the reaction rate constant kon, which 



appears in macroscopic (mass action) rate equations and which contains 
additional contributions from the diffusive behavior of the particles. For 
later reference and to make our convention for k^ clear, we give the radia- 
tion boundary condition for 3D and 2D more explicitly 

d 

47ra^gL»— G(r,i|ro,to)||r|=a„,, = Kl°Gi\r\ = a^e,t\ro,to) 



dr 
d 

27raoffi:»— G(r,i|ro,io)||r|=Qrft = K^^^G{\r\ = acs,t\ro,to). 



Clearly, (I2.14p and (I2.15P are the limit of (I2.16P when Ka -> oo and Kq -^ 0, 
respectively. The described boundary conditions render the associated solution 
nontrivial. Notwithstanding, analytical representations for the full solutions are 
known for all three cases in 3D and 2D [8 . Their numerical approximation is 
discussed in section U) In section 12.21 it will be shown how they can be used to 
remedy the issue of underestimating the number of particle encounters during a 
naive BD simulation time step. 

Other quantities of interest can be derived from the Green's function. The 
survival probability of a pair of molecules separated by Tq at time io not to react 
and thus survive until at least time t is defined by 

5'abs,rad(i|ro,io) = / Gabs.rad (r, i|ro, io)c?^r- (2.17) 

v'|r|>aotf 

Note that one has two different survival probabilities corresponding to the bound- 
ary conditions involving particle absorption. For purely reflective bcs, the sur- 
vival probability is equal to one for all times. In the cases of spherical symmetry 
that we consider here, the survival probability depends only on the radial Green's 
function g(r, i|ro,to) 

/•oo 

Sahs,rad{t\ro,to) ^ ujd gahs,iad{r,t\ro,to)r''^^^dr. (2.18) 

Thus, for all three boundary conditions the radial Green's functions can be 
obtained by 

u^dg{r, t\ro,to) = I G(r, t\ro, to)d''n (2.19) 

where d Q, means the infinitesimal surface element of the unit sphere in d di- 
mensions. 

2.2 Reversible bimolecular reactions and collision detec- 
tion via Green's functions 

In this section we introduce the key components of the proposed simulation 
algorithm. In particular, we will discuss expressions for 

• the encounter probability that will permit collision detection, 

• the propagator that will tell us how to (re-)sample the new positions of the 
molecules involved in encounters, 

• the reaction probability that allows us to determine the pairs of encountered 
molecules that subsequently undergo a chemical reaction. 



• describing unimolecular reactions, in particular back-reactions C ^^ A + B. 

To this end we again describe two-particle systems in terms of the difFusional 
behavior of a point-like particle around a sphere in 3D and 2D, but instead of 
a description based on the Fokker-Planck equation as in 12.11 we switch to a de- 
cription in terms of individual trajectories and adopt the stochastic differential 
equations' point of view. Indeed, we will emphasize the relationship between 
both approaches by recalling how the Green's functions of the initial and bound- 
ary value problem arise naturally in the trajectory picture without ever using a 
partial differential equation. Although most of the ID examples can be found in 
textbooks [19], we will present them here for two reasons. First, we think that 
the chosen representation makes the underlying physics evident and facilitates 
the derivation of the relevant expressions greatly. Second, by clarifying the role 
of the first-passage time, the relationship with event-driven approaches can be 
established, cp. \5\ 

2.2.1 Encounter Probability 

Consider a point-like particle with initial position Fq at time to in the vicinity of 
an absorbing sphere dS and let S denote the region of space internally bounded 
by dS. In a naive BD simulation, the particle's displacement r — Fq in the next 
time step is sampled according to the free-space Langevin equation (|1.2[) . i.e. 
according to a Gaussian pdf, cp. (j2.13p . 

After the simulation time step the particle may end up in the "forbidden" 
region S, which necessarily - assuming continuity of the trajectory - means that 
the particle has hit the boundary somewhere in the course of the time step. In 
the case of an absorbing boundary this is equivalent to the occurrence of an 
reaction. However, even if the particle's position after the time step is outside 
the boundary region, one cannot conclude that there was no encounter during 
At. This problem was analyzed in reference [B^ and it was shown how to address 
it by computing an " encounter probability" , given that the particle arrives at 
r after the time step and given that it started a ro. We will rederive that 
expression in a way that will be useful for later considerations. To this end we 
consider (mathematical) Brownian motion W4, which is a stochastic process, 
i.e., for every i > 0, Wt(-) : uj £ il -^ W{t,uj) is a random variable in the 
probability space {fl,A, Prob). We refrain from giving a more precise definition 
of Wi and {ft,A,Proh) that can be found elsewhere [K]. Note that at the level 
of the stochastic differential equation Wj is related to the physical stochastic 
trajectory X* by dXt = ^/l2D)dWt. 

For simplicity, let us consider a ID Brownian motion whose trajectories start 
at a point xq = at time to — without loss of generality. Furthermore, we 
assume the presence of an absorbing boundary at a; = a > xq. Therefore, the 
"allowed" region is [—00, a). The result of every simulation step corresponds to 
the stochastic outcome of an experiment and we are interested in assigning a 
probability to the following events: 

A := {there was no encounter with the boundary for all s < t}, (2.20) 

and 

B := {Xt = x\Xt„ = xo}. (2.21) 



We wish to find the conditional probability 

Prob(A|B) (2.22) 

starting from the Wiener probability distribution 

Fw{x,t\0,0) := Pi-oh{Wt < x\Wt,=o ^ 0} 

^ ^ z-y'/'^'dy (2.23) 



V2'Kt 

Henceforth, for notational simplicity, we will often suppress the condition |M^fQ=o = 
0. 

In this context, the presence of boundary conditions can be incorporated by 
introducing additional random variables. More precisely, at the level of indi- 
vidual trajectories an absorbing boundary can be realized by terminating any 
trajectory as soon as it hits the boundary for the first time. Hence, for absorbing 
boundaries the first-passage time Ta plays a central role. It is defined by [^[?T| 

Ta = mi{t e [0, t] -.Wti [-0O, a[}. (2.24) 

Note that Wr^^ = a. With the help of the first passage time one may define a 
reflected Brownian motion 

We are now in a position to calculate the joint probability distribution func- 
tion PToh{Wt < x,Ta < t). Using ([225]), {uj\Wt > 2a - x} C {uj\Ta < t}, the 
reflection principle and (|2.23p . one obtains 

ProhiWt <x,Ta<t) = Prob(M^f >2a-x;Ta< t) 
= 1 - Prob(iyt <2a-x) 

= -L= e-y'/'^dy (2.26) 

V ZTTt J2a-x 

Because of (12.231) and 

Prob(Wt < x) = Prob(H/t <x,Ta<t)+ Prob(Wt < x, r^ > t), (2.27) 

it follows 

Prob(lVt <x,Ta>t) = Prob(VFt < x) ~ Prob(VFt < x, r^ < t) (2.28) 
= -±= ( r e-y"/^'dy - r e-y"/^'dy 

The joint probability distribution function (|2.28p takes into account all Brownian 
paths that end at a; at i and whose first-passage time is bigger than i, and thus 
never encountered the boundary during the time interval. The associated pdf is 

fw{x,t\Q,Q) = ^VYOh{Wt <X,Ta> t) 
^ ¥mh{Wt =X,Ta> t) 

^ ^e-"'/2t _ ^-i.-2af/2t\ (2.29) 



/2TTt 



which is the Green's function of the initial and absorbing boundary problem of 
the ID diffusion equation. Importantly, this relation between the joint probabil- 
ity density function of the first passage time and the absorbing Green's function 
also holds in 2D and 3D. 

Because (I2.29P is Prob(A n B), we can now calculate (I2.22p . the probability 
that the Brownian mover never touched the boundary within the time step, given 
that Wt — X, according to Bayes' formula, cp. [6] 

^ , , ,^.r N Froh(Wt — x\Ta> t) 

ProHr.>tm^x) . pLb(I^,:,) ' = 

^^^f^f^^l^\ .: l-p_(x,t|.o,to) (2.30) 

(|2.30p provides a way to detect an encounter of particles, even if after the 
simulation time step the particles have no overlap. Henceforth, we refer to 
Pcnc{x,t,xo,to) as encounter probability . 

Although we considered the 3D case for simplicity, the key relationship be- 
tween first-passage time and the Green's function satisfying absorbing bcs does 
also hold in 2D and 3D and the whole line of reasoning can be extended to these 
cases. However, to this end one has to substitute in (I2.30|) the appropriate pdfs 
(Green's functions), which are known |5], cp. section HI (14.11) and (14. 4p . 

2.2.2 Propagation 

For a completely diffusion limited reaction, the encounter probability is equiv- 
alent to the reaction probability: upon detecting an encounter the algorithm 
replaces the pair of particles by its reaction products. However, if a chemical 
reaction is not completely diffusion limited, a particle reacts only with a finite 
rate upon hitting the encounter surface and is otherwise reflected. Therefore, 
the particles that did encounter each other must not be propagated according 
to the free-space Green's function, but the propagation needs to take into ac- 
count the reflective encounter, cp. [H [H US [H [111301 123 ■ In a naive BD 
simulation, reflecting bcs are frequently incorporated in an analogous manner to 
absorbing bcs: The terminal positions of only those particles that overlap after 
a time step At are reset to the boundary r G dS [53] ■ This procedure introduces 
two errors: First, as in the case of absorbing boundary conditions, the number of 
actual reflections is underestimated. Second, the resetting is only justified if the 
Brownian motion assumes a boundary position at the terminal time t — to + At, 
but never crossed the boundary during the time step. In general, however, it has 
crossed the boundary at an earlier time io < ^' < ^o + At. 

To correct for the first error, we can employ the same strategy as in the 
aborbing boundary case: We recall that (|2.30l) gives the probability that the 
particle has hit the boundary which, for an absorbing be, means that a reaction 
happened, but which means in the present context that a reflection occurred. To 
correct for the second error, instead of using the resetting method, we employ 
the Green's function satisfying reflecting boundary conditions: 

n(r) • VrG,et(r, io + Ai|ro, to)\reds = (2.31) 

However, the new positions of the encountered particles have to be sampled 
according to 

Grefl(r, to + At\ro, to) - Gabs(r, to + At|ro, to) (2.32) 

10 



instead of Gicfl(r, to + ^t|ro,io)- This can be seen as follows: Consider in ID 
the reflected Brownian motion 

^*-(r* w 'iTj-" (2-33) 

\ 2a-Wt \iWt> a ^ ' 

and the corresponding probability distribution Prob(VKt < a;) for x < a. The 
associated probability density function is the Green's function of the initial and 
reflective boundary problem of the FP equation. Now, the probability distribu- 
tion may be written as 

Prob(H^f <x)= Vmh{Wt <x,Ta<t)+ Vroh{Wt <x,Ta> t), (2.34) 

because {ra < t} and {ra > t} are mutually exclusive events and due to (|2.33p . 
As discussed in section 12.2. 1[ the second term on the right hand side yields the 
Green's function with absorbing bcs, cp. (j2.26p and (j2.28p . taking into account 
those particular trajectories of the reflected Brownian motion that did not get 
reflected. But in the simulation one only resamples the terminal points of the 
particles that did undergo an reflection. Thus, one has to sample according to 
Prob(Wt <x,Ta< t), but this is Grcf^a;, io + Ai|a;o, io) - Gabs(a;, io + At|a;o, to). 
These results can readily be extended to higher dimensions and thus we arrive 
at the following strategy to improve the accuracy of a naive BD simulation 
of reflecting particles: Sample the positions of all particles according to the 
unbounded standard process. If two particles overlap, they must have reflected 
each other and are resampled according to (J2.32I) to avoid the error made by 
resetting the particles' positions at the encounter distance. For all other pairs of 
particles we calculate the encounter probability using equation (|2.30p . If there 
was a reflection, the new positions of the particles are resampled according to 

2.2.3 Bimolecular reaction probability 

In addition, we now assume that there is a flnite probability that the particle 
has actually undergone a reaction. To take this into account the conditional 
probability that the particle has reacted upon reflection, given its initial position 
before the reflection is ro and its position after the reflection event is r, can be 
employed. Again, Bayes' formula can be used to flnd this conditional reaction 
probability, but the condition can no longer be described in terms of Gficc(r, to + 
At|ro,to) as in (J2.30L because one considers only the pairs that got reflected 
and these were propagated according to G,cfl(i", io + Ai|ro,to) — Gabs(r, io + 
At|ro,io) instead of according to the free-space Green's function. Now, it only 
remains to flnd the numerator in Bayes' formula. To this end one may argue 
in a similar way as in 12.2.21 Because only reflected pairs can react, they have 
to be described by Gi.cf(r, io + Ai|ro,io), but as pointed out in 12.2.21 the Giof 
has also contributions from the paths that did not lead to an encounter (and 
which can be subtracted by Gabs), but in the presence of partially absorbing bcs 
also contributions from paths that involved an reflection, but nevertheless did 
not react. These contributions are taken into account by Grad- Thus, for the 
numerator in Bayes' formula one may write G,cf(i", ^o + At|ro, to) — Grad(r, io + 
At\ro,to) and, in total, for the reaction probability one obtains according to 
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Bayes' formula, cp. also [23l El EH [26l [HI |30] 

Prcac(r, ro, At = „.,,.,, -T -^ — , , , . ,| -T 2.35 

Grcf(r, to + At|ro, io) " Gabs(r, io + At|ro, to) 

It is straightforward to see that the method for radiation boundary conditions 
reduces to the method for absorbing boundaries in the limit Ka ^>- oo. As 
Grad(r,io + At|ro,to) — > G„y,^(r.tn + At\rn.tn) it follows from (|2.35p that in 

this limit all particles that would have been reflected (which means that they 
reached the encounter distance) will react. This is exactly the method used in 
the case of absorbing boundary conditions. On the other hand, in the limit 
Ka — ^ one obtains Grad(r, io + At|ro, to)^'Gicf(r,^o + At\ro,to) and for the 
encounter probability (12.351) Picac(r, tq, At) -^ 0. Consequently, there will be no 
reactions and we recover the method for purely reflecting boundaries. 

2.2.4 Radial versions of Penc and Preac 

Finally, we would like to point out that instead of the full Green's functions 
including angle dependency, cp. (|4.ip . (J4.4I) . their radial counterparts (|2.19l) 
may be used for the calculation of Pe„c and Preac- At least in 3D this offers the 
advantage that one can avoid the numerical approximation of an integral, cp. 
section |4?T] In 2D the situation is different, but for small times we will derive 
an alternative expansion that does not involve an numerical integration, cp. 14.21 
Note that even if the radial versions of Ponc and Prcac are used, the particles are 
propagated according to the full propagator (|2.32l) . 

2.2.5 Unimolecular reactions 

The considerations made in the previous sections can be applied to bimolecular 
reactions. To include also the decay of single particles A — >■ reactionproducts we 
assume that this reaction can be described as a Poisson process. The probability 
that the next reaction occurs between t and t + dt is given by 

q{t\to) = fcde-^-''(*-*«)dt, (2.36) 

with kd being the particle's decay rate. To obtain the single molecule reac- 
tion probability that there was a reaction during the small but finite (i.e. not 
infinitesimal) time step At, (j2.36p has to be integrated: 

i-to+At 

Kct(At) = / q{t\to)dt = 1 - e-'^-^*. (2.37) 

J to 

Typically, one is interested in the unimolecular reaction C ^ A + B. Therefore 
the question arises how the distance between A and B after the time step should 
be chosen. To answer this we note that when C decays, the molecules A and 
B are at contact |r| — a^s and hence obviously encountered each other. We 
can thus apply the strategy adopted in the section about refiective boundary 
conditions. More precisely, we propagate the molecules according to (I2.32p with 
kol = flcff- Finally, (12.35^ is used to decide if the molecules escape or recombine. 
We conclude this section by pointing out that the methods described in this 
section may also be applied in the presence of a deterministic potential, without 
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requiring Green's function of the full Smoluchoswki equation, which can rarely 
analytically represented anyway. As described in |6|, also in this case the Green's 
functions of the pure diffusion equation can be used by separating the contribu- 
tions of the deterministic and the random force. More precisely, within one time 
step, first, the particles' displacements are calculated according to the determin- 
istic forces only and subsequently sampled according to Brownian motion. The 
encounter and reaction probabilities are then calculated by taking into account 
only starting and terminal position of the Brownian motion, instead of starting 
and terminal point of the combined motion. 

3 Simulation algorithm 

In the previous sections we have shown how Green's functions can be used to 
improve the naive BD simulation of the diffusive behavior of particles in the 
vicinity of reactive boundaries. The encounter and reaction probabilities, (j2.30p 
and (j2.35p . respectively, play a crucial role. They permit to compute the proba- 
bility that there was a reaction during a time step, even if the particles have no 
overlap after the time step. In this way, underestimating the number of actual 
reactions can be avoided. As a consequence, the convergence behavior can be 
improved and the time step can be chosen bigger than the one used in a naive 
BD simulation. The simulation algorithm consists of the following steps: 

1. Propagate all particles' positions according to the solution of the free-space 
Smoluchowski equation, i.e. for the displacements one has X^^+At — ^to = 
v2DAtN{0,l), where N(0,1) denotes the unit normal random variable 
with vanishing mean and variance equal one. 

2. Check for all overlapping particles if the involved species pairs imply ab- 
sorbing, reflecting or radiation boundary conditions. 

(a) In case of absorbing bcs, the particles underwent a reaction and are 
replaced by the reaction products. 

(b) In case of reflecting boundaries, the particles' position is reset by the 
following two steps. First, sample the R coordinate (j2.2|) according to 
(|2.13p . Second, the relative distance vector (|2.3p is sampled according 
to Grcfi(r, to + Ai|ro, to) - Gabs(r, to + At|ro, to) <^M- 

(c) In case of radiation boundary conditions, proceed in the same way as 
in the case of reflecting boundary conditions. In addition, calculate 
(|2.35p . where r is given by the new, according to (|2.32p resampled, 
positions. Draw a uniform random number ^. If ^ < Proac, there was 
a reaction. Replace the pair of particles by the corresponding reaction 
products. 

3. Calculate for all pairs of particles that do not overlap the encounter prob- 
ability (|2.30p . Draw a uniform random number ^. If ^ < Penc, there was 
an encounter. Apply to all in this way found encounter pairs the steps 2 
(a), (b), (c). 

4. Calculate for each molecule that can undergo a unimolecular reaction the 
associated survival probabilty (I2.37p . Draw a uniform random number ^. 
If ^ < Prcacj tli6 molecule decayed and is replaced by its reaction products. 
Propagate the reaction products and check for escape or recombination as 
described in section 12.2.21 12.2.31 
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5. Increase the system time by At. Repeat steps 1., 2., 3, 4. 

Details 

• In step (3a), a particular molecule might possess several neighbors which 
fulfill the encounter criterion ^ < Ponc- In this way the algorithm generates 
a list of encounter candidates for a particular molecule. To decide which 
of the candidate molecules has actually been involved in the encounter 
process one can proceed in a way similar to how the Gillespie algorithm 
decides which reaction channel is going to "fire" [16]: More precisely, if 
there are A'cand neighbors which fulfill the encounter criteria, one can apply 
the following algorithm 

— Calculate the sum of the (already in step 3 calculated survival proba- 
bilities Pane) of all candidate molecules. 

(3.1) 









Ptot 


= 


A'cand 


Pcnc,i 




— Calculate the 


ratios 


Ei 


Ptot 


J 


= 1, 


• ■ • -^*cand- 


It follows 


o< ^^ 


-1 Pcnc, 


u 


eU 


7^cnc,i 




nd 

i^enc,i 



(3.2) 

Ptot Ptot Ptot 

— Sample a uniform random number ^. Find the j which satisfies: 

Z^i Penc.i , t ^ Z^j Penc.i /„ „n 

Ptot Ptot 

— Pick the j + 1 molecule as " encounter" molecule and continue with the 
algorithm starting from 3 (a). 

Propagation by G'rcfl(r,io + At|ro,to) - Gabs(i",io + At|ro,to) means the 
following: 

— First, the radial propagator 

^dr'^^^ {gicfi{r, to + At\ro, io) - 5abs(?', io + Ai|ro, to)), 
which gives the probability of finding the radial coordinate in the inter- 
val [r, r + dr[ at time to + At, is used to sample a new radial coordinate. 

— Then, ^(Grefl(r, 8, to + At|ro, Oo, to) - Gabs(r, 8, to + At|ro, Go, to)) , 
which yields the probability of finding the angle between r and ro 
inside the interval [0,0-1- dQ[, given that the radial coordinate is r 
at t, is used to draw a new angle. The factor J' is defined as J^ := 
27rsin(0)r^ in 3D and as J^ := r in 2D, respectively. 

To decrease computational cost one can introduce a cut-off radius rcut-off 
such that only neighboring particles within that radius are tested for en- 
counters according to described procedure. Following [B] the cut-off radius 
might be determined by requiring that the survival probability (j2.17p sat- 
isfies 

S{to + At|reut-off , to) = 10-3. (3.4) 

for the chosen simulation time step At. Depending on the type of the be 
that models the interaction between the considered species pair, one uses 
either Sabs or S'rad- 
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Finally, we would like to emphasize again that the described algorithm is 
not exact, but relies on approximations, as is obvious in the case of several 
encounter neighbors described above. The accuracy will be analyzed in detail in 
a forthcoming publication. 

4 Numerical approximation of the Green's func- 
tions 

4.1 General aspects 

The crucial ingredient of the suggested algorithm is the use of the Green's func- 
tions, which permit to increase the simulation time step. On the other hand, 
the calculation of the Green's functions is the most costly part of the algorithm. 
Hence, an efficient calculation of the Green's functions is required. 

The analytical representation for the Green's function describing the diffusion 
of a point particle around a partially absorbing sphere is known [Q . In a scaled 
form, which is suitable for numerical approximations, it is given in 3D by 

oo 

e^'"''' F^+i/2{R,x)F^+i/2{Rq,x)x dx. (4.1) 
The functions F^j, where v = n + 1/2, are defined by 

F (R t) = (2h''^ + l)[J^(Rx)Y,Ax)-Y^iRx).J,Ax)]-2xlJARx)Y;{x)-YJRx)jUx)] . . „^ 
"^ ' '^ {[{2h3-D + i)J^(x)-2x,J'^{x)]^ + [{2h3'0 + l)Y^{x)-2xYl{x)]^}^/^ ' ^ ' 

Here, R = r/a, Rq = r^j a denote the dimensionless relative radial coordinates 
after and before the time step, respectively, and denotes the angle between 
the corresponding relative position vectors. Furthermore, t = Dt/o? is the 
dimensionless time and 

/l3D := /i3D„ .^ J^ (4,3) 

might be thought of as a dimensionless reaction constant. J„+i/2,i^n+i/2 ^'^^ 
the fractional Bessel functions of first and second kind [T] , respectively, and P„ 
denote the Legendre ploynomials of order n [T]. For h^^ —> cx) and h^^ —>■ 0, (j4.ll) 
reduces to the Green's functions satisfying absorbing and reflecting boundary 
conditions, respectively. 

In 2D the corresponding Green's function is 

G,^d{R:0,t\Ro) = —^ V cos(n0) / e-^-''Cn{R,x)Cn{Ro,x)xdx (4.4) 
where the functions Cn{R, x) are defined by 

^ MRx)[xy;{x) - fe^py„(x)] - y^{Rx)[xj:,{x) - l^^Ux)] 

' ([xJaa;)-/i2i?j„(a:)]2 + [a;y^(a;)-^2Dy^(a;)]2)i/2 ' ^'^ 
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and 



h^"" := h 



2D, 



,.2D 



27rD' 



(4.6) 



Note that in (j4.4p the Bessel functions of first and second kind are of integer order 
[T]. Again, one can obtain the Green's function for absorbing and reflecting bcs 
from (|4.4p by taking h"^^ -^ cx) and ft^° -^ 0. 

The 3D and 2D radial Green's functions are given, up to a factor, by the ze- 
roth term of the expansions (|4.ip and (j4.4p . respectively. Now, a special property 
of fractional Bessel functions is that they can be expressed in terms of elementary 
functions [I] eq. 10.1.11, 10.1.12] 



Jl/2ix) 
Yl/2ix) 



sin(a::) 



cos(a;) 



(4.7) 
(4.8) 



Thus, the integral that involves the function F1/2 and that yields the radial 
Green's functions can actually be solved explicitly for all three boundary condi- 
tions [S]. Here we give them again a scaled form 



and 



^(exp 



{R~Ror 

4r 



exp 



{R^+Ro-2y 
4r 



(4.9) 



5rad(-R, T|i?o) 



8Tia'->RRo ' 



exp 



(-B.--Ro)^ 
4t 



exp 



{R+Ro-2) 



/Attt exp (k^t + {R + Rq - 2)k) erfc (k^/t - 



iR+Ro-2) 
2^7 



For convenience, we put 



ifP + inaD 



(4.10) 



(4.11) 



47raD 
gicf is given by the same expression as ^rad with k^^ — 0. 

By contrast, in the 2D case the corresponding integrals involve Bessel func- 
tions of integer order for which no expression in terms of elementary functions 
is known and hence, the calculation of the 2D radial Green's functions requires 
a numerical integration. As a consequence, the 2D simulation is more costly in 
this regard than a 3D simulation, because, as described in section |3l the radial 
Green's functions are used by a random number generator to sample new po- 
sitions. Moreover, the radial versions of Pcnc and Prcac mentioned in 12.2.41 that 
permit faster simulations in 3D, loose this advantage in the 2D case. In par- 
ticular for smaller r this leads to a substantially increased computational cost 
when compared to the 3D case. To address this problem we will derive in section 
4.21 a small time expansion for the radial Green's functions satisfying absorbing, 
reflecting and radiation boundary conditions, respectively. 

The analytical representations (14.11) and (j4.4p can be numerically approxi- 
mated by the use of three-term recurrence relations [lOj , which are satisfied by 
the functions P„, cos(n6'), J^, Y^: For P„ one has [II eq. 8.5.3] 



(n + l)P„+i{x) = (2n + l)xPn{x) - nP„-i(x), 



(4.12) 
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where Po{x) — l,Pi{x) — x. Both the Bessel functions Ji,{x),Y^{x) satisfy the 
same three-term recurrence relations, but due to the fact that for Ju{x) one has 
to iterate down it is more convenient to write them as [1] eq. 9.1.27] 

Mx) = ^(^^±llj,+i(a;) - J,+2(.t) (4.13) 

x 

Y,+2{x) = ^^^^illy,+i(a;) - Y,{x) (4.14) 

X 

Note that (|4.13p can be used for both Bessel functions of integer and fractional 
order, cp [B]. Finally, the c„ := cos{n9) satisfy 

Cn+2 = 2 cos(a;)c„+i - c„, (4.15) 

where cq = l,ci = cos(a::). For the integration procedure we chose a Gaussian 
integration rule. The advantage of this method is that not only the weights and 
abscissas {xi}, but also the values of the Bessel functions at the abscissas (but 
not at Rxi, R^Xi) can be calculated in advance, i.e. before the actual simulation 
starts. 

4.2 Small time expansion of 2D Green's functions 

For small times, i.e. in terms of the dimensionless time r := Dt/a^ < 0.01, the 
representation of the Green's functions given above becomes cumbersome. The 
involved integrals become increasingly difficult to solve numerically due to the 
oscillatory character of the Bessel functions which become less and less damp- 
ened by the exponential factor. Furthermore, as already discussed above, in 2D 
the integrals yielding the radial Green's functions cannot be solved analytically, 
in contrast to the 3D case. As a consequence, the 2D case is actually more cum- 
bersome than the 3D case in the context of the suggested simulation algorithm 
(actually in the context of any algorithm based on the 2D Green's functions). 

To address this problem we will show in the following how one can obtain an 
alternative expression for the 2D radial Greens's functions in terms of a small 
time expansion and thus, effectively, solve the integral. 

[S] demonstrates how the Laplace transform technique can be used to find 
small time expansions for a number of solutions to the heat equation. To the best 
of our knowlege, the specific small time expansions we are deriving here have not 
been discussed yet. We make the following ansatz for the Laplace transform of 
the radial Green's function that satisfy certain boundary conditions [8] 

9{r, q\ro) = gircc{r, g|ro) + ghc{r, q\ro) (4.16) 

Here 

^-'=('^''^l^")-2^1/o(<zr)i^o(<Z.o) r<r, ^^.17) 

is the Laplace transform of the radial free-space Green's function, q is defined by 
q := ^/%, where p denotes the Laplace domain variable. The part ^bc that takes 
into account the boundary condition is a solution to the Laplace transformed 2D 
diffusion equation 

-T^-H 1 q ghc = 0- (4.18) 

dr^ r ar 
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The general solution to (|4.18p is AKo{qr) + Blo(qr), where Io{x), Ko{x) refer 
to the modified Bessel functions of first and second kind, respectively, and of 
order zero ]1\. Because we require lima;_>.oo ffbc — >■ 0, and lini2;_j.oo -^o (2;) — >■ 00, 
B has to vanish and hence, 5bc('', <z|^o) = AKo{qr). A is determined by the 
requirement that the complete Green's function g{r,q\ro) (|4.16p satisfies either 
absorbing, reflective or radiation boundary conditions, respectively 

9{r,q\ro)\r=a = (4.19) 

dg{r,q\ro) 



dr 



n=a = (4.20) 



dg{r,q\ro) ^^ 

— \r=a = h g{r,q\ro)\r=a- (4.21) 

h^^ has been defined in (I4.6p . In the following we will use h instead of h^^ . 
We switch to dimensionless variables q :— qa, h :— ha, R := j'/a, Rq := ro/a 
and obtain from (HTTC)) . (|IT7)) . (HTTg|) . (|I^ and (|i?^ for the corresponding 
ghciR,(i\Ra) components 

5abs(i?,9|i?o) = --^^^Ko{qRo)KoiqR). (4.22) 

~grei{R,q\Ro) - "^:§^^'o(«~^o)i^o(9i?). (4.23) 

1 hUq)-I'(q) 
5rad(i?,<?|i?o) - - . „,-; ^i^o(«i?o)i^o(9i?). (4.24) 

2^^|i^o(9)-i^^(g) 

Instead of applying the inversion theorem for the Laplace transformation to the 
complete solution of the boundary problem given above (|4.16p . we make the 
following detour. First, we note that in almost all for the algorithm required 
expressions the free-space part cancels out, so we focus on the boundary com- 
ponents. Second, for small times we are interested in obtaining expansions in 
powers of p~^, or equivalently, q~^. To this end we exploit the asymptotic ex- 
pansions of the modified Bessel functions for large arguments [T] 



Mkiv) 




(4.25) 
(4.26) 



fc=0 



where the coefficients are 



(4^^-P)(4^^-3^)...(4r.^-(2fc-l)^) 
ak[i^) = ^k (4.27) 

Using these asymptotic expansions and 

i;,ix) = h{x), (4.28) 

K'oix) = -K^{x) (4.29) 

one can transform the exact expressions (|4.22l) . (|4.23p and (|4.24p . according to 
the rules for multiplication and division of asymptotic series [7, , to the following 



18 



asymptotic expansion in powers of g ^ . 

^^^^ -^ "i^^^vW — ^ — ^„ 9- ' ^'-''^ 

^-^'-^ -^ i^WM^ ^ to ^^ ■ 

Note that the different boundary conditions lead to the same generic form of their 
asymptotic expansions, which only differ with regard to sign and the expansion 
coefficients 7fe. For convenience, we give the first two expansion coefficients 

abs ^ 2i?i?o - (i? + Rq) 



SRRq 
^b. _ 9(i?2 + i?2) + 2RRo - 4R Ro{R + Ro)+AR^R'^ 





aos _ ^ ■ u^ ' " "' 

^2 f28i?2i?2 



,ad _ i?i?o(6 + f6/t) + i? + i?o 

8-R-Ro 
,^d _ 9(i?^ + i^g) + 2RRo + i?i?o(i? + i?o)(12 + 32/i) + i^^j^a^gg ^ ^Q2h + 256fe^ 
^2 ~ 128i?2i?2 

71°^! 72^^^ can be obtained from 7^^'*, 72^'^ for h = 0. Wc sec that the explicit form 
of the coeffcients becomes quickly cumbersome, so that in a simulation all the 
required coefficients are calculated by iterative use of (|4.27l) . 

To find the small time expansion in the time domain we use the Laplace 
transforms [S] 

g-qx fD^^''^ 



-^ ( — ) e-3ro (4.32) 



a; 



4-— -> (4i)«/2i"erfc(^^) (4.33) 

pl+n/2 ^ ' ^2VDi 

The functions i"erfc(a;) are defined by 

/•OO 

i"erfc(a;) ;= / i''-^eric{^)d^ (4.34) 

J X 

and 

i°erfc(a:) :=erfc(x), (4.35) 

cp.[H]. Moreover, 

1 1 _ 2 

i erfc(x) := ierfc(x) := — ^e ^ — a;erfc(a;). (4.36) 

None of the integrals in (j4.34p have to be calculated, because the i"erfc(x) func- 
tions satisfy the recursion relation 

2n i"erfc(x) = i''-2erfc(x) - 2x i"-ierfc(a;) (4.37) 
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that allows to calculate the (J4.34I) swiftly. In this way, we finally arrive at the 
expressions for small times in the time domain 



5abs(-R, r|i?o) ~ - 
and 

5rof, rad(^,'''|-Ro) "" 



1 



^■Ka'^y/TRR^ 



TT-i/^e- 



{R+B.Q-2y 



-. OO 



7r(fi,^o)i"-'erfc 



R + Ra-2 



i^v^y 



(4.38) 



1 



Ana^y/^RRo 



1 /o (R+rtp-^r 



iE7r-^^(i?,i?o)i"-^erfc(^±^ 

n— 1 ^ ^ 



i^v^y 



(4.39) 



Again, the differences between the expressions for the different boundary condi- 
tions manifest only in the sign and the 7„ coefficients. 

Let us reconsider the reaction probability Preac expressed in terms of the 
radial 2D Green's function to relate it with the radiation boundary constant. 
Using the derived short time expansions and [5] 



/^e-\Mx) = - - ^ 



(4.40) 



1-h 



2t 



(4.41) 



one obtains for sufficiently small times 

gradiR, t\Ro) - gabs(-R, t\Ro) _ 
greiiR,T\Ro)-gMR,r\Ro) ~^ '^\R + Ro-2 

and hence, the reaction probability defined in section [2.2.31 can be expressed to 
first order in the time step as 

r 2r t 



Pr. 



R + Ro-2 



' Tra{r + tq — 2a) 



(4.42) 



Thus, for sufficiently small times the reaction probability is the ratio of the 
radiation boundary reaction constant Ha to the ratio of the area a{r + tq — 2a) 
and the time step. 

Finally, we would like to make three comments. First, for large h, the quotient 
hK° (a)-'' K'l'^a) ^^ ^^^ accurately represented by its asymptotic expansion. In 
these cases it is advantageous to start from 



1 



h lo(qa) 

h Ko(qa) 



We can expand (|4.43p as a series in h 



q 



3.ad - ^K,{qr,)K,[qr)^l{l ^ 



hiqa) Ki{qa) 



loiqa) Koiqa) 



q\^ Ki{qa) 
Ko{qa) 



(I) 



Ia{qa) K^iqa) 



K 



M\ 



2 r 



Ko{qa)J 



h{qa) Ki(qa) 
h{qa) Koiqa) 



(4.43) 



(4.44) 

+ ...}. 

(4.45) 
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Still, there can be ranges of h where both expansions fail. In general, however, 
these ranges are small and do not pose a serious obstacle to the described method. 

Second, although we only considered the gbc parts, small time expansions for 
the complete Green's functions can easily be found by using the known expres- 
sions for the free-space Green's functions, cp (j4.16p and (|4.17p . 

Third, in principle, the small time expansion technique described here for the 
integrals involving the function C„=o (|4.5p can also be applied for the integrals 
with C„^o- However, for larger n the convergence behaviors of asymptotic ex- 
pansions worsen, especially for r, tq close to a. Still, these expansions can prove 
useful for very small time steps and r, tq 3> a. 

5 Appendix: Relation to first-passage time sim- 
ulation algorithms 

As mentioned in the introduction, the presented algorithm may be regarded as a 
coarse-grained version of an event-driven first-passage time simulation algorithm. 
Here, we are demonstrating this relationship explicitly. The line of reasoning in 
the sections 12.2.11 \TIT^ and [^.2.31 emphasizes the crucial role of the first-passage 
time. Neglecting events {Wt < x,Ta < t} leads to the \/At error in the naive BD 
simulation. The conditional probability penc can be used to detect these events, 
so that the remaining error in the formalism given above is due to the uncertainty 
with regard to the first-passage time. The detection of an encounter event via 
Pcnc solely provides an upper bound, i.e. Axa — At. On the other hand, the 
exact knowledge of the first-passage time would allow to set up an event-driven 
algorithm, with time steps equal to sampled first-passage time. It turns out that 
one can determine the first-passage time with any desired accuracy and without 
explicitly using first-passage time distributions, but only with the expressions 
derived so far. A similar construction has been described in [29] to determine 
the last reflection time. 

Consider a ID Brownian motion, describing a molecule in the vicinity of an 
boundary with Wtg = 0. After one simulation time step it is found at Xf = 
Wtg+At- As we have just seen in (J2.30I) . Ponc can be employed to decide if there 
was at least one encounter. Assuming now that there was at least one encounter, 
it follows to < ''a < ^0 + At. The idea is now by combining iterative bisection 
of the time interval [to, to + At] and application of Pcnc to approximate the 
first-passage time Ta with in principle arbitrary precision. More precisely, one 
considers the time intervals [to,to + At/2] and [to + At/2,to + At]. Application of 
(j2.30p requires the reconstruction of the intermediate point x„i = Wtg+At/2 the 
molecule assumed at t ^ to + At/2, given that initially it was at Wtg — Xi and 
ended up at Wt^+At — Xf. Obviously, it is not correct to sample the intermediate 
position according to the free-space Green's function. Instead, one considers the 
increments Ci = Wto+At/2 - Wt^, Cm = Wto+At - Wta+At/2 and C/ = d + Cm 
which are by definition Gaussian random variables with (d) = (^m) = and 
(^2^ = (^^) = a^ := At/2, hence (C/) = and {^j) = 2a'^ = At. Thus, the 
conditional probability density is again a Gaussian where the first moment equals 
half the distance between initial and terminal point 

P^^' - ^™l?/ - ^f> - p(^^ = ^^) - ^^ exp( ^2 ) 15.1) 
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By sampling according to (j5.ip (or the higher dimensional equivalent) one can 
construct the intermediate point. Once the intermediate point is known, one 
can use (|2.30p to decide if there was an encounter in [to, to + At/2]. If there was 
an encounter, it follows that io < '''a < ^o + At/2 and we repeat the described 
procedure for the intervals [to, At/2] and [xi,Xrn]- If there was no encounter, 
it follows that to + At/2 < Tq < to + At and we repeat the procedure for the 
intervals [At/2, At] and [x„i,x/]. 

Thus, the chosen time step of the simulation corresponds to a cut-off of the 
first-passage time uncertainty and, conceptually, one might consider a family 
of the presented algorithm {SimAlgJAt, "indexed" by different time steps as a 
sequence of renormalization group transformations. In this picture, the limit- 
ing case with vanishing uncertainty corresponds to an exact event-driven first- 
passage time formalism. We would like to point out that these considerations 
might not only be of pure conceptual interest. It is conceivable that the given 
method could be used to resolve certain time periods of the simulation with finer 
detail. Moreover, using the algorithm as part of a hybrid algorithm, simulations 
on finer and coarser scales could be done in a more controlled way. 
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